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1. Introduction 

inge of applications for random lattices spam 

J 



There is a large range of applications for random lattices spanning from materi- 
als science,EJ modelling of the large-scale structure of the universe,0 kinetic growth 
models,cl sand pileso to quantum field theoryE A random lattice is a way of dis- 
cretizing a system without introducing any kind of anisotropy. It can, however, 
also be used to describe cellular structures, such as grain mosaics, biological tissues 
and foamsfl'Qtm So it is important to have efficient algorithms for generating ran- 
dom lattices, and, furthermore, algorithms which allow the random lattice to be 
dynamically maintained throughout a simulation. 

An extensively used type of random lattice is the Delaunay random lattice.0 
which is a triangulation of space based on a given set of N points. The Voronoi 
tessellation^ is dual to this lattice and gives a division of space into cells. It is de- 
scribed as follows: For each point one determines the region of space which is closer 
to this point than to any other point. This procedure called the Voronoi tessellation 
of space will divide the space into convex cells. There exist other ways of divid- 
ing space into cells, such as the Laguerre partitionjlj and, recently a generalized 
Voronoi construction was used in order to model soap froths.! 
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In the present paper we discuss static and dynamic algorithms for generating 
Delaunay random lattices determined according to the Voronoi tessellation of space. 
We will restrict the discussion of the Voronoi construction to two dimensions; it 
can be generalized to any dimension. We determine the computation time for 
our algorithms, and will denote this by T(N) for a random lattice with N points. 
Results for both the Poissonian random lattice (PRL)E3 and the vectorizable random 
lattice (VRL)H will be given. The PRL is obtained when the points are put at 
independently and randomly chosen positions all over the entire lattice and then 
connected according to the Voronoi construction. This implies that the PRL random 
lattice is completely free of anisotropics]! 

The VRL lattice is a modified form of the PRL random lattice. Now the space 
is divided into square cells and each cell contains exactly one point. The position of 
a point inside a cell is chosen randomly. The motivation for introducing the VRL 
random lattice was to have a vectorizable program when solving an equation on 
the lattice, e.g., the Laplace equatiorLdescribing the numerical growth of Diffusion- 
Limited Aggregation (DLA) clusters.B Furthermore, the VRL lattices possess for all 
practical purposes no anisotropics and since a vectorized code can increase compu- 
tation time significantly VRL lattices are very useful. 

First, in Sec. 2, the Voronoi and Delaunay constructions are described in more 
detail. In Sec. 3 we discuss our different algorithms. Section 4 contains a discussion 
of the performances and finally we conclude. 

2. Voronoi Tessellations and Delaunay Random Lattices 

The Voronoi construction or tessellation for a given set of points is defined as follows: 
For all points determine the associated cell consisting of the region of space nearer 
to this point than to any other point. Whenever two cells share an edge they are 
considered as neighbours. By drawing a link between the two points associated 
to (located in) the cells one obtains the triangulation of space that is called the 
Delaunay lattice. The Delaunay random lattice is dual to the Voronoi tessellation 
in the sense that points correspond to cells, links to edges and triangles to the 
vertices of the Voronoi tessellation. 

The triangulation (tessellation) of the plane with Np points (cells), Nl links 
(edges) and Nt triangles (vertices) are constrained by the Euler relation 

Np — Nl + Nt — X (1) 

where x is the Euler characteristic, which equals 2 for a graph on a surface with 
the topology of a sphere and for a torus. Furthermore, 3Nt = 2iV^ holds for a 
triangulation (a triangle has three sides/links and 3iVr then counts each link twice). 

An equivalent way to obtain the Delaunay random lattice is the following: Con- 
struct for all sets of three points the circle defined by these points. If there are no 
other points inside this circle the three points will define a Delaunay triangle and 
are connected by three links. The centers of these circumscribing circles will be the 
vertices of the Voronoi tessellation.!] 
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3. Algorithms 

Three different algorithms for generating two-dimensional Delaunay random lattices 
and the associated Voronoi tessellations are presented. The programs were written 
in the C programming language and run on workstations. 

First we describe an algorithm for the VRL case, and, secondly, the PRL random 
lattice. Thirdly, a dynamic algorithm for generating both the VRL and PRL lattices 
is described. We call this algorithm dynamic because it allows one to change a 
random lattice slightly and then obtain the new random lattice by only carrying 
out local rearrangements of the links. Examples of VRL and PRL random lattices 
are shown in Fig. 0. 

3.1. Vectorizable random lattices 

In the case of the vectorizable random lattice the region where the points are dis- 
tributed is first divided into square cells as shown in Fig. |l](a). Each cell is of unit 
area. Then exactly one point is put into each cell with the position inside the cell 
being random and uniformly distributed. 

Each point only can be connected to points in the 36-cell neighbourhood as 
shown in Fig. |l|(a). In Fig. |l](b) the probabilities for the center point to be con- 
nected to another point in this neighbourhood are shown (cf. Ref. 3). So if a point 
belongs to a triangle the other two points are contained in two of the cells of this 
neighbourhood. Therefore, to find a triangle one only has to search for neighbour- 
ing points in this restricted area. An upper limit for the distance two neighbouring 
points can have is \/4 2 + 2 2 = \/20 ~ 4.472 as can be seen from Fig. p](a). 



Fig. 1. (a) The division into cells in order to obtain the VRL random lattice. Also shown is the 
initial triangle ABC together with the 36-cell neighbourhood of the initial point A (see text). For 
the front link AB a point P is searched such that the circle through these three points does not 
contain any other point. The triangle APB is then a Delaunay triangle, (b) Probabilities for a 
point to be connected to another point in the 36-cell neighbourhood for a point. 

3.1.1. VRL algorithm 

After putting the random points the algorithm determines an initial triangle (cf. 
Fig. 0(a)). One point A is chosen. For each pair of points in the neighbourhood of 
A it is checked if the circle defined by these three points contains any other point 
of the neighbourhood. The first triplet of points found for which this circle does 
not contain any other points is connected with three links. These links are given a 
counter clockwise orientation, as shown by the arrows in Fig. [j](a). 

When the lattice is completed each link is either connected to two triangles or 
it is a border link, i.e., it is part of the border which limits the entire lattice. Each 
link is assigned a flag front which has the value if the two triangles connected 
with this link have been found or if the link is a border link. If only one triangle has 
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been found front has the value 1. We call such a link a front link, due to the fact 
that this is the front of the search algorithm. So after finding the initial triangle 
its three links AB, BC and CA are front links. The following steps are now carried 
out: 

For each front link AB check if a point P in the neighbourhood of A can be 
found such that the circumscribing circle of the triangle APB does not contain any 
other point (cf. Fig. 0(a)) . This is then a Delaunay triangle. First it is checked if 
the point P is on the right or the left side of the link AB, i.e., if the determinant of 
the matrix with the column vectors AB and AP is negative or positive. a Then the 
next point in the neighbourhood has to be checked. If no point P completing AB 
to a Delaunay triangle can be found the link AB is a border link and its front flag 
is assigned the value 0. The program continues by treating the next front link. If a 
point P with this property is found there are three new links AP, PB and BA. One 
of them BA is identical with the just checked front link AB, it only has a different 
orientation. The front flag of AB is assigned the value 0. The probability to find 
a connected point to P is not equal for all cells of the neighbourhood but is much 
smaller for the outer cells. We use this by beginning the search at the cells with the 
highest probability (cf. Fig. |l|(b)) and then going outwards to the cells with lower 
probability of connection. 

It is checked if the new links AP and PB have already been found. If not, they 
are new front links. Otherwise they already were saved as front links and their 
front flag is assigned the value 0. One notices that it is not necessary to look 
through the entire array of links to check whether AP and PB have already been 
found. In our algorithm we use an array connect in which we save for each point 
the links connected to this point. Therefore it is only necessary to check whether 
the links AP and PB are already saved in the connect arrays of the points A and B. 
If there is still a front link left treat the next one, otherwise the Delaunay random 
lattice is finished. 

Finally the Voronoi tessellation is determined. For each point P all triangles 
which have P as one site are ordered counterclockwise. The centers of the circum- 
scribing circles of the triangles are connected. The polygon constructed in this way 
is the Voronoi cell which is dual to the point P. In summarized form the algorithm 
looks like the following: 

• Put random points into the region 

• Construct the initial triangle 

• Construct the Delaunay random lattice. Do the following steps as long as 
there are front links AB left: 

1. Find point P which completes front link AB to a Delaunay triangle. If 
not possible, AB is a border link. 

a If P is on the left side, APB cannot be a new triangle of the random lattice because in this case 
it has already been found. 
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2. Check if the two new links are already stored as front links. If not they 
are new front links. 

• Construct the Voronoi tessellation by connecting the centers of the circum- 
scribing circles of all triangles connected to a point 

The computation time T(N) for constructing the VRL random lattice is pro- 
portional to the number of points TV. This is due to the fact that for each point the 
neighbour points can only be found in the 36-cell neighbourhood. 

3.2. Poissonian random lattices 

In this section we describe an algorithm for generating a random lattice consisting 
of points with uniformly distributed x and y values. In the limit of infinite lattices 
this corresponds to a Poisson process (see, e.g., Ref. 13). 

The actual generation of a random lattice based on a Poissonian distribution of 
lattice points has many aspects in common with the approach used in the previous 
Section 3.1. Different techniques are, however, necessary in some places mainly 
due to the fact that for the PRL the number of neighbours is not bounded. The 
algorithm used here is described in detail in Ref. 4 and based on the one proposed 
in Ref. 5. 

3.2.1. PRL algorithm 

During the initialization process, a starting triangle is identified and its edges are 
defined to be active links, similar to the front links of the previous section. There- 
after the list of active links is treated until there are no more active links left. To 
identify an active link a flag treated is used, which is initialized to zero for every 
new link found. To generate the random lattice one has then to go through the 
following steps: 

• Take the next active link: 

1. Draw a circle around this link and check if it contains points 

2. If so, choose the one which leads to the largest radius 

3. If not, move the center of the circle perpendicular to the link by a small 
amount e dependent on the number of points present in the region. Draw 
a new circle. Repeat until the circle contains points. Select the one which 
leads to the smallest radius 

• Store the point found as a new site and create two new links by connecting it 
to the two sites of the considered link. 

• Check if these new links already exist. If they do, remove the active attribute, 
if not append them as new active links 

• Remove the active attribute from the considered active link 
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Fig. 2. Finding a new point to construct a new triangle starting from one link of an existing 
triangle. Mi , M 2 and M 3 are centers of the circles, the C's are points of the lattice. 



One has to search for every active link the next neighbour in the sense of the 
Delaunay construction. A way to do this is shown in Fig. |^. A circle is drawn 
through the endpoints of the link AB, its center being the link's midpoint Mi. One 
has then to check, if there are points inside that circle (e.g., C\). If there are, one 
has to choose the point which gives the largest radius by drawing a circle determined 
by the point and the endpoints of the link. If in the beginning there were no points 
in the circle, one has to move the center of the circle by a small amount e dependent 
on the density of points into the direction where the neighbour should be located 
(e.g., M 2 , M 3 ), until one or more points are inside the circle (e.g., C2, C3). One 
chooses then the point giving raise to the circle with the smallest radius. This 
method is used because testing if a point is inside a circle or not is much faster than 
the radius calculation. 



3.2.2. Data structure 

The two basic elements of the random lattice, sites and links, are organized as a 
linked list of special data structures which allow to store the information needed: 



SITE{ 



POINT p; 


/* 


Coordinates */ 




int border; 


/* 


Border Site Flag */ 




int nNb; 


/* 


Number of Neighbours 


*/ 


NBNQDE *pFirstNbNode ; 


/* 


Address of the First 


Neighbour Site */ 


SITE *pNextSite; 


/* 


The Next Site in the 


List */ 


SITE *pNextS!Site; 


/* 


The Next Site in the 


Current Bin */ 



LINK{ 



POINT pi, p2; 


/* 


Coordinates of Constituting Points */ 


SITE *pSitel,*pSite2; 


/* 


Addresses of these Points */ 


int orient ; 


/* 


Orientation of the Link */ 


char border; 


/* 


Border Link Flag */ 


char treated; 


/* 


Treated Link Flag */ 


LINK *pNextLink; 


/* 


The Next Link in the List */ 



}; 



The site list is initialized only once in the beginning. During the construction of the 
lattice neighbours are, however, dynamically appended to each site. The address of 
the first neighbour is given by pFirstNbNode. Each appended neighbour is again 
a structure containing information about the address of the neighbour site and the 
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link connecting it: 
NBN0DE{ 

NBNODE *pNextNbNode ; /* The Next Neighbour Site */ 

SITE *pNb; /* Address of the Current Neighbour */ 

LINK *pLink; /* Address of the Link Connecting */ 

>; /* the Site with this Neighbour */ 

In contrast to the arguments given in Section 3.1 one does not know where ex- 
actly to find the neighbour site. For the search of this neighbour it is however not 
acceptable to do a complete search of all points, since this would imply a time be- 
haviour proportional to TV 2 . Alternatively one could use a tree-like structure, where 
one can identify points in the neighbourhood. Another possibility is to implement 
binning, i.e., divide the lattice into sublattices so that in each bin one finds on the 
average the same number of sites. This method has been implemented here. The 
pointer pNextSlSite contains the address of the next site in a considered bin. Thus 
for each link — on the average — only a fixed number of sites has to be searched. This 
number should be equal to 5-10 sites per bin, while the 8 neighbouring bins arc 
searched for each active link. 

During the construction of the lattice the algorithm will encounter situations, 
where a neighbour found gives raise to new links which have already been found, 
and are still waiting in the active link list to be treated. In this case one has to 
mark these links as treated by setting the corresponding flag treated in the link 
structure. Treated links are skipped when processing the link list. 

The stopping criterion depends on the problem one is investigating. If there are 
no periodic boundary conditions, one has to introduce a limit beyond which the 
algorithm should stop to look for new points. If the problem has periodic boundary 
conditions, no special stopping criterion is needed. 

3.2.3. Speed and memory 

The speed of a program depends on the setup of the algorithm and the realization 
on a given computer. The best performance one could have is 0(N) 1 i.e., the order 
of the computation time T(N) increases linearly with the number of points treated. 

To avoid any supcrflous calculation the following measures have been taken in 
this algorithm: 

• While creating the lattice one has to check whenever a new neighbour and 
therefore two links are found, that these links are not yet in the list of active 
links. To minimize the impact of this task on the time behaviour, each site 
knows through the structure NBNODE not only about its neighbours pNb but 
also about the connecting links pLink. So only the links already attached to 
the site have to be searched, a number which is constant on the average. 

• Searching potential neighbours, one has to draw circles and move their centers. 
Care is taken, that the quantity e by which the centers are moved depends 
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on the density of points. If one would not account for that, more and more 
heavy radius calculations would be necessary as the density grows. 

With these measures taken the algorithm should have a time behaviour like O(N), 
since for every point — on the average-a fixed number of calculations and memory 
allocations are carried out. 

These theoretical considerations do, however, not reflect reality. Dynamic mem- 
ory allocation is intuitively and easily programmable; the computer will on the other 
hand be obliged to maintain longer and longer lists of allocated memory blocks, so 
memory allocation time will grow faster than the system size. 

That is why in this program memory is allocated in huge blocks from the start. 
Whenever new memory is needed (for new neighbours or links) it is taken from 
these "pools". The computer thus has no additional work to do for maintaining 
memory allocation tables. b Allocating memory in the beginning is in fact possible 
since one knows that the average number of neighbours for a site is 6. Particularly 
for large lattices (N > 10,000) this is very advantageous, since memory management 
tends to become a substantial part of the calculation time. Allocating memory from 
the start one has the other advantage that the linking of the lists is not necessary 
anymore, and this saves memory. 

3.3. Dynamic Random Lattices 

Now we describe a method to calculate a random lattice using a Dynamic Random 
Lattice algorithm (DRL), which makes it possible to dynamically maintain the 
random lattice in a simulation.EH The idea is to start from an existing lattice (e.g., 
a triangular lattice) and then successively change this lattice by moving the points 
in order to obtain the final random lattice. The lattice fulfils at all times the Voronoi 
construction and this ensures that also the final lattice will be a Delaunay random 
lattice. 



Fig. 3. Random lattices, (a) and (b) VRL lattices, (c) and (d) PRL lattices, with periodic boundary 
conditions obtained with the DRL algorithm. The lattices contain 400 points and were obtained 
by starting from a triangular lattice. 

By the DRL algorithm, the VRL random lattice is obtained by moving the 
points — one by one — to new positions randomly chosen in the cells (cf. Fig. |l|(a)) 
associated to the points (see Figs. |](a)-(b)). The PRL random lattice is obtained 
by moving the points to new positions randomly chosen over the entire lattice (see 
Figs. S(c)-(d)). The time dependence for the DRL algorithm is more complicated 
than for the other two algorithms due to the "search" part of the algorithm (see 
below) and will have the slowest performance in the case of the PRL random lattice 
and the fastest in the case of the VRL random lattice. 

b The total time for computing large lattices will increase more than linearly since the computer 
will be forced to swap memory pages. 
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3.3.1. DRL algorithm 

We now describe the basic elements of the algorithm (see Refs. 7 and 14 for a 
more detailed description). When a point is moved from one place to another it 
is possible to view all the changes of the lattice as link-flip processes. A link-flip 
process is when, e.g., the link EP in Fig. ^(a) is flipped to DF as in Fig. [|(b). 
In this case it is a neighbour-loosing process for the points P and E, whereas it is 
a neighbour-gaining process for the points D and F. Intuitively, when a point is 
moved, one follows the trajectory of the point and tests whether a link-flip process 
will take place. The condition for this to occur is: i) either the point P enters 
the circumscribing circle of a triangle (corresponding to the points A, B and C 
in Fig. |J(a)), which then implies a neighbour-gaining process (cf. Fig. |J(b)) or, 
ii) the point P leaves the circumscribing circle for three points (corresponding to 
the points D, E and F in Fig. ||(a), which do not define a triangle), and this implies 
a neighbour-loosing process (cf. Fig. [|(b)). 

Fig. 4. Part of a Delaunay random lattice (a) before and (b) after the point P has been moved. 
The link-flip process appearing to the left of the point P is a neighbour-loosing process whereas 
the process to the right is a neighbour-gaining process. The thick lines represent the links of the 
random lattice whereas the thin lines are the edges of the Voronoi cells. 

Even though it is impossible to predict the number of processes that will take 
place, the full set of changes can always be viewed as consisting of the neighbour- 
gaining and loosing processes taking place by following the point along its trajectory. 
It is, however, not necessary to carry out the updating of all the processes along 
the points trajectory. When a point has to be moved one first removes it from the 
lattice leaving a lattice with N — 1 points, and updates the neighbourhood. Then 
the point is put at its new position, and the location of this new position, i.e., to 
which triangle it belongs in the lattice, is determined. Since one does not know into 
which of the existing triangles on the random lattice the point is put, the way this 
triangle is determined is by a search through the lattice (see Ref. 14). Finally, the 
point is added and the neighbourhood is updated by carrying out neighbour-gaining 
and loosing processes. In schematized form the DRL algorithm looks like this: 

• Generate a triangular lattice with N points and the topology of a torus (pe- 
riodic boundary conditions) 

• For all points P on the lattice do: 

1. remove the point P and the links to its neighbours 

2. calculate/update the new links 

3. put the point P at a new (randomly chosen) position 

4. search through the Delaunay random lattice with N—l points and locate 
in which triangle the new point belongs 
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5. add the point P and calculate/update the new links followed by neigh- 
bour flipping processes (cf. Fig. 0) 



3.3.2. Data structure 

The random lattice algorithm is implemented using a dynamic pointer structure. 
For each point a SITE structure (see below) is defined where information associated 
to the site is stored. Included in this information is a nearest-neighbour pointer list 
NN. For each site in this list there is a pointer referring to the location of the original 
site itself, so that the information associated to a site is only stored once. The sites 
and nearest neighbour linked lists are organized as follows: 



SITE { 








POINT 


p; /* 


coordinates */ 




int 


n; /* 


point number */ 




int 


q; /* 


coordination number */ 




NN 

}; 


*pNN; /* 


pointer to first Nearest Neighbour 


*/ 


NN { 








SITE 


*pSite; 


/* pointer to the site */ 




CIRCLE circle; 


/* circumscribing circle */ 




int 


pos ; 


/* position relative to the site 


*/ 


NN 


*pNN; 


/* pointer to the Next Neighbour 


*/ 



}; 

The updating of the lattice when a flip process takes place is carried out by 
rearranging the pNN pointers in the nearest-neighbour pointer lists. When a new 
triangle appears or an old one disappears information on the circumscribing circles 
is updated and stored in the NN nearest- neighbour pointer list (cf. Ref. 14). 

3.3.3. Performance 

The average computation time of the DRL algorithm can be written as 

T(N) = 0(N{R{N) + S(N) + A(N))) (2) 

where R(N) is the average time to remove a point, S(N) the search time and A(N) 
the time to add a point. For a triangular lattice the coordination number is 6 so 
R(N) and A(N) would both imply 6 nearest-neighbour operations. From Eq. ([!]) 
follows that the average coordination number for a triangulation is 6, so both R(N) 
and A(N) will on average be approximately constant. 

The search time SSJSf) is proportional to the distance between the old and new 
position of the point. EJ In the VRL case the new positions of the points are within 
a fixed distance from the old positions, i.e., S(N) is a constant, and we expect the 
total time T(N) to scale linearily with N, i.e., an O(N) behaviour. For the PRL 
lattice the search will in the worst case be proportional to the size of the system, 
i.e., \/N, leading to an 0(N>/N) behaviour for the computation time. 
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4. Results 

All simulations we report in this section were carried out using Sun Sparc 10 work- 
stations. In Table [l] and Fig. |^ are shown the computation times for the static VRL 
and PRL algorithms as well as for the dynamic algorithm, denoted here as DVRL 
and DPRL. For the static version one notices an approximate linear increase in 
computation time as a function of the number of points N on the lattice. For large 
systems the performance is worse due to the fact that almost all the memory of the 
computer is used. In the VRL case the expression T(N) w 0.207V milliseconds fits 
the data very well whereas the PRL algorithm is a little slower with the expression 
T(N) w 0.247V milliseconds providing a good fit. c 

Table 1. The upper part shows the performance of the VRL algorithm for the VRL random lattice 
and for the PRL algorithm the PRL random lattice. In the lower part the values for the dynamic 
algorithm are given for the VRL and PRL random lattice, denoted here as DVRL and DPRL. The 
times are measured in seconds. 



# Points, N 


100 


1000 


10,000 


50,000 


100,000 


200,000 


300,000 


VRL 


0.04 


0.18 


1.72 


8.90 


18.3 


38.6 


60.3 


PRL 


0.02 


0.20 


2.05 


11.3 


23.0 


47.3 


73.9 


DVRL 


0.11 


0.99 


10.15 


50.6 


101.0 


207 


308 


DPRL 


0.14 


1.70 


29.20 


264.8 


713.8 


1938 


3476 



Fig. 5. Computation time for the VRL, PRL and DRL algorithms in order to generate random 
lattices with N points. The times are measured in seconds. Also shown are curves representing a 
linear increase and an Nt/N increase in computation time. 

For the DRL algorithm one notices a linear relation between execution time and 
the number of points TV for the VRL random lattice, T(N) w l.OliV milliseconds. 
In the PRL case the behaviour asymptotically approaches the 0(Ny/~N) behaviour, 
as discussed in Sec. 3.3.3. d An expression of the form T(N) = aN(b + y/N) (with 
a ps 150 ms and b w 65) follows our curve (cf. Fig. ||), and this has the anticipated 
0(Ny/~N) behaviour for large N. When in advance one has a knowledge of the 
distribution of the points it may be possible to reduce the search time by using 
"bins" (cf. Sec. 3.2.2), but we have only investigated the general DRL algorithm 
without using any information about the location of the points. 

5. Conclusions 

In the present paper we have discussed different algorithms for generating Delaunay 
random lattices and the associated (dual) Voronoi tessellations. We have described 

c The linear behaviour of the algorithm is only valid if every bin contains about the same number of 
points. For distributions of x and y values which are not uniform, the binning should be adjusted. 
" The initialization time for generating the initial triangular lattice where memory for the dynamic 
pointer structure is allocated is included in the above computation times. This initialization time 
scales approximately linearily as cN with c 0.157 ms. 
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how the algorithms have been implemented in the C programming language using 
dynamic pointer structures. 

For the vectorizable random lattice and the Poissonian random lattices we found 
that the computation time is roughly linear. When the system size approaches the 
memory limit of the computer the time dependence increases more rapidly. We 
also discussed a dynamic random lattice algorithm. This algorithm allows for a 
random lattice to be maintained during a simulation and can also be used in order 
to obtain a random lattice but with a slower time performance than the VRL and 
PRL algorithms. 
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